package commons.utils;

import commons.utils.primitive_geom.Point_d;
import commons.utils.primitive_geom.Rect_d;
import commons.utils.primitive_geom.Segment_d;
import commons.utils.primitive_geom.Vector_d;
import gnu.trove.list.array.TIntArrayList;
import java.util.ArrayList;
import java.util.Collections;



public class UtilsGeom {

	public static Point_d getCentreGravite(ArrayList<? extends Point_d> vPoint) {
		if (vPoint.isEmpty())
			return null;
		
		Point_d res = new Point_d(0, 0);
		for (Point_d pt : vPoint) {
			res.x += pt.x;
			res.y += pt.y;
		}
		res.x /= vPoint.size();
		res.y /= vPoint.size();
		return res;
	}
	
	public static Point_d getClosestPoint (Point_d pt, ArrayList<Point_d> vPoints) {
		double closestDist = Double.POSITIVE_INFINITY;
		Point_d closestPoint = null;
		int nbPointToTest = vPoints.size();
		for (int i = 0; i <nbPointToTest ; i++) {
			Point_d pt_toTest = vPoints.get(i);
			double dist = pt.distTo(pt_toTest);
			if (dist<closestDist) {
				closestDist = dist;
				closestPoint = pt_toTest;
			}
		}
		return closestPoint;
	}

	public static Double getDistancePtSegment (Point_d ptA, Point_d ptB, Point_d ptI) {

		// Attention au cas où le point ne se projette pas sur le segment [A,B] !
		//        -> ->
		// (1) Si AI.AB <= 0 alors la distance la plus courte est ||AI||
		//              -> ->
		// (2) Sinon si BI.BA <= 0 alors la distance la plus courte est ||BI||
		// (3) Sinon, I se projète perpendiculairement sur le segment [A,B], et la distance correspondante est calculée

		// Cas (1)
		double dxAI = ptI.x-ptA.x;
		double dyAI = ptI.y-ptA.y;
		double dxAB = ptB.x-ptA.x;
		double dyAB = ptB.y-ptA.y;
		if (dxAI*dxAB + dyAI*dyAB <= 0)
			//return Math.sqrt ( (ptA.x-ptI.x)*(ptA.x-ptI.x) + (ptA.y-ptI.y)*(ptA.y-ptI.y) );
			return Math.sqrt ( dxAI*dxAI + dyAI*dyAI );

		// Cas (2)
		double dxBI = ptI.x-ptB.x;
		double dyBI = ptI.y-ptB.y;
		double dxBA = ptA.x-ptB.x;
		double dyBA = ptA.y-ptB.y;
		if (dxBI*dxBA + dyBI*dyBA <= 0)
			//return Math.sqrt ( (ptB.x-ptI.x)*(ptB.x-ptI.x) + (ptB.y-ptI.y)*(ptB.y-ptI.y) );
			return Math.sqrt ( dxBI*dxBI + dyBI*dyBI );

		// Cas (3)
		//Area = |(1/2)(x1y2 + x2y3 + x3y1 - x2y1 - x3y2 - x1y3)|   *Area of triangle
		//Base = V((x1-x2)²+(x1-x2)²)                               *Base of Triangle*
		//Area = .5*Base*H                                          *Solve for height
		//Height = Area/.5/Base

		// !! Cette formule est instable peut provoque un phenomene de cancellation (cf. http://www.irisa.fr/sage/jocelyne/cours/precision/insa-1200.pdf page 30)
		//double area = Math.abs (.5 * (ptA.x * ptB.y + ptB.x * ptI.y + ptI.x * ptA.y - ptB.x * ptA.y - ptI.x * ptB.y - ptA.x * ptI.y));
		
		// Cette formule est stable :
		double double_area = Math.abs( dxAB*dyAI - dyAB*dxAI );
		
		double bottom = Math.sqrt ( dxAB*dxAB + dyAB*dyAB );
		return double_area / bottom; // = height
	}
	
	// idem en plus rapide si le carré de la distance suffit (pour un interclassement ou un tri de distances, par ex.)
	public static Double getDistancePtSegment_sq (Point_d ptA, Point_d ptB, Point_d ptI) {
		double dxAI = ptI.x-ptA.x;
		double dyAI = ptI.y-ptA.y;
		double dxAB = ptB.x-ptA.x;
		double dyAB = ptB.y-ptA.y;
		if (dxAI*dxAB + dyAI*dyAB < 0)
			return dxAI*dxAI + dyAI*dyAI;

		double dxBI = ptI.x-ptB.x;
		double dyBI = ptI.y-ptB.y;
		double dxBA = ptA.x-ptB.x;
		double dyBA = ptA.y-ptB.y;
		if (dxBI*dxBA + dyBI*dyBA < 0)
			return dxBI*dxBI + dyBI*dyBI;

		// !! Cette formule est instable peut provoque un phenomene de cancellation (cf. http://www.irisa.fr/sage/jocelyne/cours/precision/insa-1200.pdf page 30)
		//double area = Math.abs (.5 * (ptA.x * ptB.y + ptB.x * ptI.y + ptI.x * ptA.y - ptB.x * ptA.y - ptI.x * ptB.y - ptA.x * ptI.y));
		
		// Cette formule est stable :
		double double_area = Math.abs( dxAB*dyAI - dyAB*dxAI );
		
		double bottom_sq =  dxAB*dxAB + dyAB*dyAB;
		return (double_area*double_area) / bottom_sq;
	}

	
	public static Point_d getClosestPt (Segment_d seg, Point_d ptI) {
		return getClosestPt(seg.pt1, seg.pt2, ptI);
	}
	

//	// Renvoie le point le plus proche de I parmi le segment [A,B]
//	// Renvoie nulle si le point est POSSIBLEMENT plus loin du seg que distMaxSearch
//	public static Point_d getClosestPt (Point_d ptA, Point_d ptB, Point_d ptI, double distMaxSearch) {
//		if (ptI.x-ptA.x > distMaxSearch || ptI.x-ptB.x > distMaxSearch ||
//			ptA.x-ptI.x > distMaxSearch || ptA.x-ptI.x > distMaxSearch ||
//			ptI.y-ptA.y > distMaxSearch || ptI.y-ptB.y > distMaxSearch ||
//			ptA.y-ptI.y > distMaxSearch || ptA.y-ptI.y > distMaxSearch )
//			return null;
//		return getClosestPt(ptA, ptB, ptI);
//	}

	// Renvoie le point le plus proche de I parmi le segment [A,B]
	public static Point_d getClosestPt (Point_d ptA, Point_d ptB, Point_d ptI) {

		// Attention au cas où le point ne se projette pas sur le segment [A,B] !
		//        -> ->
		// (1) Si AI.AB < 0 alors le point le plus proche est A
		//              -> ->
		// (2) Sinon si BI.BA < 0 alors la distance la plus courte est ||BI||
		// (3) Sinon, I se projète perpendiculairement sur le segment [A,B], et la distance correspondante est calculée

		// Cas (1)
		double dxAI = ptI.x-ptA.x;
		double dyAI = ptI.y-ptA.y;
		double dxAB = ptB.x-ptA.x;
		double dyAB = ptB.y-ptA.y;
		if (dxAI*dxAB + dyAI*dyAB <= 0) { // correction juillet 2013 (avant c'était '<')
			//System.out.println("cas 1");
			return ptA;
		}

		// Cas (2)
		double dxBI = ptI.x-ptB.x;
		double dyBI = ptI.y-ptB.y;
//		double dxBA = ptA.x-ptB.x; //TODO_OK PERF : = -dxAB -> OK
//		double dyBA = ptA.y-ptB.y; // idem -> OK
		if (-dxBI*dxAB - dyBI*dyAB <= 0) { // correction juillet 2013 (avant c'était '>')
			//System.out.println("cas 2");
			return ptB;
		}


		// Cas (3)
		// Soit J le point de projection
		// J est tel que :
		// -> ->          -> ->   ->
		// JI.AB = 0  et  AI^AB = 0
		double xJ, yJ;
		if (ptA.y == ptB.y) {
			xJ = ptI.x;
			yJ = ptA.y;
		} else {
			// TODO_OK PERF reprendre les dxA ... -> OK
			// TODO_KO PERF precalculer delta dans chaque segment [A,B]
		    double Delta = (dxAB) / (dyAB);
			yJ = (Delta*Delta*ptA.y + Delta*(dxAI) + ptI.y) / (Delta*Delta+1);
			xJ = (yJ-ptA.y) * Delta + ptA.x;
		}
		//System.out.println("cas 3");
		return new Point_d(xJ, yJ);
	}
	
	
//	// Renvoie dans tabResult, le point de seg1, puis le point de seg2 tels que ces 2 points forment
//	//  la plus petite distance entre les 2 segments
//	public static double getClosestPt_seg_seg (Segment_d seg1, Segment_d seg2, Point_d[] tabResult) {
//		
//		// Cas special ou les 2 segments s'intersectent
//		Point_d ptIntersect = intersectionSegments(seg1, seg2);
//		
//		if (ptIntersect != null) {
//			tabResult[0] = ptIntersect;
//			tabResult[1] = ptIntersect;
//			return 0;
//		}
//		
//		// Sinon, on fait toutes les combinaisons de distance d'une extrémité d'un seg à l'autre seg
//		//  et on choisit la plus petite distance
//		Point_d pt1, pt2;
//		double dist;
//		
//		Point_d pt_1_1 = getClosestPt(seg1, seg2.pt1);
//		double dist_1_1 = pt_1_1.distTo(seg2.pt1);
//		
//		pt1 = pt_1_1;
//		pt2 = seg2.pt1;
//		dist = dist_1_1;
//		
//		
//		Point_d pt_1_2 = getClosestPt(seg1.pt1, seg1.pt2, seg2.pt2, dist);
//		if (pt_1_2!=null) {
//			double dist_1_2 = pt_1_2.distTo(seg2.pt2);
//
//			if (dist_1_2 < dist) {
//				pt1 = pt_1_2;
//				pt2 = seg2.pt2;
//				dist = dist_1_2;
//			}
//		}
//			
//		
//		Point_d pt_2_1 = getClosestPt(seg2.pt1, seg2.pt2, seg1.pt1, dist);
//		if (pt_2_1!=null) {
//			double dist_2_1 = pt_2_1.distTo(seg1.pt1);
//
//			if (dist_2_1 < dist) {
//				pt1 = seg1.pt1;
//				pt2 = pt_2_1;
//				dist = dist_2_1;
//			}
//		}
//				
//				
//		Point_d pt_2_2 = getClosestPt(seg2.pt1, seg2.pt2, seg1.pt2, dist);
//		if (pt_2_2!=null) {
//			double dist_2_2 = pt_2_2.distTo(seg1.pt2);
//
//			if (dist_2_2 < dist) {
//				pt1 = seg1.pt2;
//				pt2 = pt_2_2;
//				dist = dist_2_2;
//			}
//		}
//		
//		tabResult[0] = pt1;
//		tabResult[1] = pt2;
//		return dist;
//	}
	
	
	// Les segments ne doivent pas s'intersecter ni être superposés, mais peuvent avoir une extrémité commune
	//  si extrémité commune et bNoDist0=true, la plus petite distance non nulle est renvoyée
	public static double getClosestPt_seg_seg (Segment_d seg1, Segment_d seg2, Point_d[] tabResult, boolean bNoDist0) {
		assert(areSegmentInternallyIntersected(seg1, seg2) == false);
		assert(areSegmentSuperposed(seg1, seg2) == false);
		
		// 1 extrémité en commun ?
		Point_d ptCommun = seg1.getCommonExtremity(seg2);
		if (ptCommun != null) {
			if (!bNoDist0) {
				// On autorise des distance nulle
				tabResult[0] = ptCommun;
				tabResult[1] = ptCommun;
				return 0;
			} else {
				// On n'autorise pas de distance nulle
				// Recherche des extrémités libres
				Point_d extrLibre_seg1, extrLibre_seg2;
				if (seg1.pt1.equalsXY(ptCommun))
					extrLibre_seg1 = seg1.pt2;
				else
					extrLibre_seg1 = seg1.pt1;
				if (seg2.pt1.equalsXY(ptCommun))
					extrLibre_seg2 = seg2.pt2;
				else
					extrLibre_seg2 = seg2.pt1;
				// Puis on fait les 2 combinaisons entre une extrémité d'un seg et l'autre segment, on prendra alors la dist min
				Point_d pt1 = UtilsGeom.getClosestPt(seg1, extrLibre_seg2);
				double dist_seg1_extrSeg2_sq = pt1.distTo_sq(extrLibre_seg2);
				
				Point_d pt2 = UtilsGeom.getClosestPt(seg2, extrLibre_seg1);
				double dist_seg2_extrSeg1_sq = pt2.distTo_sq(extrLibre_seg1);
				
				if (dist_seg1_extrSeg2_sq < dist_seg2_extrSeg1_sq) {
					tabResult[0] = pt1;
					tabResult[1] = extrLibre_seg2;
					return Math.sqrt(dist_seg1_extrSeg2_sq);
				} else {
					tabResult[0] = extrLibre_seg1;
					tabResult[1] = pt2;
					return Math.sqrt(dist_seg2_extrSeg1_sq);
				}
			}
		}
		
		// Cas standard : pas d'extrémité en commun
		// On fait toutes les combinaisons de distance d'une extrémité d'un seg à l'autre seg
		//  et on choisit la plus petite distance
		// on ne teste plus bNoDist0, car on est sur que la distance ne sera pas nulle
		Point_d pt1, pt2;
		double dist_sq;
		Point_d pt_1_1 = getClosestPt(seg1, seg2.pt1);
		double dist_1_1_sq = pt_1_1.distTo_sq(seg2.pt1);
		pt1 = pt_1_1;
		pt2 = seg2.pt1;
		dist_sq = dist_1_1_sq;
		
		Point_d pt_1_2 = getClosestPt(seg1, seg2.pt2);
		double dist_1_2_sq = pt_1_2.distTo_sq(seg2.pt2);
		if (dist_1_2_sq < dist_sq) {
			pt1 = pt_1_2;
			pt2 = seg2.pt2;
			dist_sq = dist_1_2_sq;
		}
			
		Point_d pt_2_1 = getClosestPt(seg2, seg1.pt1);
		double dist_2_1_sq = pt_2_1.distTo_sq(seg1.pt1);
		if (dist_2_1_sq < dist_sq) {
			pt1 = seg1.pt1;
			pt2 = pt_2_1;
			dist_sq = dist_2_1_sq;
		}
				
		Point_d pt_2_2 = getClosestPt(seg2, seg1.pt2);
		double dist_2_2_sq = pt_2_2.distTo_sq(seg1.pt2);
		if (dist_2_2_sq < dist_sq) {
			pt1 = seg1.pt2;
			pt2 = pt_2_2;
			dist_sq = dist_2_2_sq;
		}
		
		tabResult[0] = pt1;
		tabResult[1] = pt2;
		return Math.sqrt(dist_sq);
	}


	//**************************************************************************
	//**************************************************************************
	//**************************************************************************

	public static double getAbsCosinusBetweenSeg (Segment_d seg1, Segment_d seg2) {
		double produitLongueur = seg1.getLength() * seg2.getLength();
		if (produitLongueur==0)
			return 1;
		double cos = ((seg1.pt2.x-seg1.pt1.x) * (seg2.pt2.x-seg2.pt1.x) + (seg1.pt2.y-seg1.pt1.y) * (seg2.pt2.y-seg2.pt1.y)) / (produitLongueur);
//		if (cos<-1 || cos>1)
//			cos = cos+1-1;
		assert (cos>=-1.00001 && cos <=1.00001) : cos;
		return Math.abs(cos);
	}
	
	// Les segments doivent partager un point commun qui est ptCommun
	// Cela permet de calculer l'angle sans avoir à orienter correctement les segments
	// Angle compris dans [0 ; 2.Pi[
	public static double getAngle_betweenSegment(Point_d ptCommun, Segment_d segA, Segment_d segB) {
		assert (segA.pt1.equalsXY(ptCommun) || segA.pt2.equalsXY(ptCommun)) && (segB.pt1.equalsXY(ptCommun) || segB.pt2.equalsXY(ptCommun));
		
		double xCenter = ptCommun.x;
		double yCenter = ptCommun.y;

		double xA = segA.pt1.x!=xCenter ? segA.pt1.x : segA.pt2.x;
		double yA = segA.pt1.y!=yCenter ? segA.pt1.y : segA.pt2.y;

		double xB = segB.pt1.x!=xCenter ? segB.pt1.x : segB.pt2.x;
		double yB = segB.pt1.y!=yCenter ? segB.pt1.y : segB.pt2.y;

		double lgA = Math.sqrt((xA-xCenter)*(xA-xCenter)+(yA-yCenter)*(yA-yCenter));
		double lgB = Math.sqrt((xB-xCenter)*(xB-xCenter)+(yB-yCenter)*(yB-yCenter));
		if (lgA==0 || lgB==0)
			return 0;

		double scalaire = (xA-xCenter)*(xB-xCenter) + (yA-yCenter)*(yB-yCenter);
		double z_vectoriel = (xA-xCenter)*(yB-yCenter) - (yA-yCenter)*(xB-xCenter);

		double angle;
		if (z_vectoriel>=0)
			// si le vecteurB est dans la moitié supérieur du cercle trigo basé sur le vecteurA
			//  alors acos donne la bonne valeur de l'angle
			angle = Math.acos(scalaire/(lgA*lgB));
		else
			// sinon, acos dans l'angle du vecteur symétrique de B par rapport à l'axe horizontal
			angle = 2*Math.PI - Math.acos(scalaire/(lgA*lgB));

		return angle;
	}

	// Détermine si les 4 segments pris 2 à 2 (0+1 et 2+3) se croisent réellement en JCenter et ne constituent
	//  pas une sorte de "double V tête-bêche"
	public static boolean isThereRealCroisement(Point_d jCenter, Segment_d seg0, Segment_d seg1, Segment_d seg2, Segment_d seg3) {
		// Le test est : par rapport à seg0, le vecteur qui doit faire le plus petit angle et celui qui doit faire le plus grand angle
		//   sont seg2 et seg3 ou l'inverse (mais pas seg1)
		double angle01 = getAngle_betweenSegment(jCenter, seg0, seg1);
		double angle02 = getAngle_betweenSegment(jCenter, seg0, seg2);
		double angle03 = getAngle_betweenSegment(jCenter, seg0, seg3);

		return (angle02<=angle01 && angle01<= angle03) || (angle03<=angle01 && angle01<= angle02);
	}



	// Determine si les rect englobant des segment s'intersectent
	public static boolean areSegPotentiallyIntersected(Point_d seg1_pt1, Point_d seg1_pt2, Point_d seg2_pt1, Point_d seg2_pt2) {
		
		if (seg1_pt1.x < seg1_pt2.x) {
			if (seg2_pt1.x < seg2_pt2.x) {
				if (seg1_pt2.x<seg2_pt1.x || seg1_pt1.x>seg2_pt2.x) return false;
			} else {
				if (seg1_pt2.x<seg2_pt2.x || seg1_pt1.x>seg2_pt1.x) return false;
			}
		} else {
			if (seg2_pt1.x < seg2_pt2.x) {
				if (seg1_pt1.x<seg2_pt1.x || seg1_pt2.x>seg2_pt2.x) return false;
			} else {
				if (seg1_pt1.x<seg2_pt2.x || seg1_pt2.x>seg2_pt1.x) return false;
			}
		}
		if (seg1_pt1.y < seg1_pt2.y) {
			if (seg2_pt1.y < seg2_pt2.y) {
				if (seg1_pt2.y<seg2_pt1.y || seg1_pt1.y>seg2_pt2.y) return false;
			} else {
				if (seg1_pt2.y<seg2_pt2.y || seg1_pt1.y>seg2_pt1.y) return false;
			}
		} else {
			if (seg2_pt1.y < seg2_pt2.y) {
				if (seg1_pt1.y<seg2_pt1.y || seg1_pt2.y>seg2_pt2.y) return false;
			} else {
				if (seg1_pt1.y<seg2_pt2.y || seg1_pt2.y>seg2_pt1.y) return false;
			}
		}
		return true;
	}


	// Renvoie true si les segments seg1 et seg2 sont intersectes (y compris s'ils ont une extrémité en commun)
	public static boolean areSegmentIntersected (Segment_d seg1, Segment_d seg2) {
		return intersectionSegments(seg1, seg2) != null;
	}
	public static boolean areSegmentIntersected (Point_d seg1_pt1, Point_d seg1_pt2, Point_d seg2_pt1, Point_d seg2_pt2) {
		return intersectionSegments(seg1_pt1.x, seg1_pt1.y, seg1_pt2.x, seg1_pt2.y, seg2_pt1.x, seg2_pt1.y, seg2_pt2.x, seg2_pt2.y) != null;
	}
	
	// Superposition complete ou partielle
	// ce concept est tres liee a la precision des calculs aussi on considere que les segments sont superposes si c'est le cas en dessous
	//   d'un certain seuil
	private static double TOLERANCE_SUPERPOSITION = 0.00001;
	private static double TOLERANCE_SUPERPOSITION_sq = TOLERANCE_SUPERPOSITION*TOLERANCE_SUPERPOSITION;
	public static boolean areSegmentSuperposed (Segment_d seg1, Segment_d seg2) {
		return areSegmentSuperposed(seg1.pt1, seg1.pt2, seg2.pt1, seg2.pt2);
	}
	public static boolean areSegmentSuperposed (Point_d seg1_pt1, Point_d seg1_pt2, Point_d seg2_pt1, Point_d seg2_pt2) {
		// la regle est simple, il faut que parmi les 4 extremites des segments, au moins 2 appartiennent
		//  a l'autre segment s'ils n'ont aucune extremite en commun
		//  3, s'ils ont 1 extremite en commun, et 4 s'ils ont 2 extremite en commun (mais le resultat est trivial alors)
		
		int nbExtremitEnCommun = 0;
		if (seg1_pt1.distTo_sq(seg2_pt1) < TOLERANCE_SUPERPOSITION_sq) nbExtremitEnCommun++;
		if (seg1_pt1.distTo_sq(seg2_pt2) < TOLERANCE_SUPERPOSITION_sq) nbExtremitEnCommun++;
		if (seg1_pt2.distTo_sq(seg2_pt1) < TOLERANCE_SUPERPOSITION_sq) nbExtremitEnCommun++;
		if (seg1_pt2.distTo_sq(seg2_pt2) < TOLERANCE_SUPERPOSITION_sq) nbExtremitEnCommun++;
		
		int nbPtOnSeg = 0;
		int nbPtOnSeg_expected;
		if (nbExtremitEnCommun==0)
			nbPtOnSeg_expected = 2;
		else if (nbExtremitEnCommun==1)
			nbPtOnSeg_expected = 3;
		else
			return true;
		
		
		if (getDistancePtSegment_sq(seg2_pt1, seg2_pt2, seg1_pt1) < TOLERANCE_SUPERPOSITION_sq)
			nbPtOnSeg++;
		if (getDistancePtSegment_sq(seg2_pt1, seg2_pt2, seg1_pt2) < TOLERANCE_SUPERPOSITION_sq)
			nbPtOnSeg++;
		if (nbPtOnSeg>=nbPtOnSeg_expected)
			return true;
		if (getDistancePtSegment_sq(seg1_pt1, seg1_pt2, seg2_pt1) < TOLERANCE_SUPERPOSITION_sq)
			nbPtOnSeg++;
		if (nbPtOnSeg>=nbPtOnSeg_expected)
			return true;
		if (getDistancePtSegment_sq(seg1_pt1, seg1_pt2, seg2_pt2) < TOLERANCE_SUPERPOSITION_sq)
			nbPtOnSeg++;
		if (nbPtOnSeg>=nbPtOnSeg_expected)
			return true;
		return false;
	}

	// Renvoie true si les segments seg1 et seg2 sont rellement intersectes et pas simplement connectés entre eux
	public static boolean areSegmentInternallyIntersected (Segment_d seg1, Segment_d seg2) {
		if (seg1.hasCommonExtremity(seg2)) return false;
		return intersectionSegments(seg1, seg2) != null;
	}

	public static boolean areSegmentInternallyIntersected (Point_d seg1_pt1, Point_d seg1_pt2, Point_d seg2_pt1, Point_d seg2_pt2) {
		if (seg1_pt1.equalsXY(seg2_pt1) || seg1_pt1.equalsXY(seg2_pt2) || seg1_pt2.equalsXY(seg2_pt1) || seg1_pt2.equalsXY(seg2_pt2)) return false;
		return intersectionSegments(seg1_pt1.x, seg1_pt1.y, seg1_pt2.x, seg1_pt2.y, seg2_pt1.x, seg2_pt1.y, seg2_pt2.x, seg2_pt2.y) != null;
	}


	/**
	 * Renvoie l'intersection des segments 1 et 2, null si pas d'intersection ou infinité
	 */
	
	// renvoie null si l'intersection est sur les extrémité
	public static Point_d intersectionSegments_interne(Segment_d seg1, Segment_d seg2) {
	    if (seg1.hasCommonExtremity(seg2)) return null;
		return intersectionSegments (seg1.pt1.x, seg1.pt1.y, seg1.pt2.x, seg1.pt2.y, seg2.pt1.x, seg2.pt1.y, seg2.pt2.x, seg2.pt2.y);
	}
	public static Point_d intersectionSegments(Segment_d seg1, Segment_d seg2) {
	    return intersectionSegments (seg1.pt1.x, seg1.pt1.y, seg1.pt2.x, seg1.pt2.y, seg2.pt1.x, seg2.pt1.y, seg2.pt2.x, seg2.pt2.y);
	}
	/**
	 * Renvoie l'intersection des segments [AB] et [CD], null si pas d'intersection ou infinité
	 */
	public static Point_d intersectionSegments(Point_d ptA, Point_d ptB, Point_d ptC, Point_d ptD) {
	    return intersectionSegments (ptA.x, ptA.y, ptB.x, ptB.y, ptC.x, ptC.y, ptD.x, ptD.y);
	}

	// on arrondi a 5 decimale le resultats, car des fois les erreurs de calcul
	// peuvent donner des choses bizares
	// si les segment sont superposés (meme partiellement) et qu'ils ont une extrémité en commun, alors
	//  l'une de ces extrémités en commun sera retournée (-> il y a un autre test si on veut reellement tester la superposition)
	public static Point_d intersectionSegments(double xA, double yA, double xB, double yB, double xC, double yC, double xD, double yD) {

		//*** Test trivial des extrémités communes *****************************
		if (xA==xC && yA==yC)
			return new Point_d(xA, yA);
		if (xA==xD && yA==yD)
			return new Point_d(xA, yA);
		if (xB==xC && yB==yC)
			return new Point_d(xB, yB);
		if (xB==xD && yB==yD)
			return new Point_d(xB, yB);
		
		//*** Test trivial des rectangles englobants ***************************

		//if (Math.max(xA, xB) < Math.min(xC, xD)) return null;
		//if (Math.min(xA, xB) > Math.max(xC, xD)) return null;
		//if (Math.max(yA, yB) < Math.min(yC, yD)) return null;
		//if (Math.min(yA, yB) > Math.max(yC, yD)) return null;

		// equivalent de ce qui precede, mais avec 6 à 8 tests seulement au lieu de 12
		// et moins d'appels de fonctions
// PLUS LENT QUE CE QUI SUIT
//		double xMin_AB, xMax_AB, yMin_AB, yMax_AB;
//		double xMin_CD, xMax_CD, yMin_CD, yMax_CD;
//		if (xA < xB) {
//			xMin_AB = xA;
//			xMax_AB = xB;
//		} else {
//			xMin_AB = xB;
//			xMax_AB = xA;
//		}
//		if (yA < yB) {
//			yMin_AB = yA;
//			yMax_AB = yB;
//		} else {
//			yMin_AB = yB;
//			yMax_AB = yA;
//		}
//		
//		if (xC < xD) {
//			xMin_CD = xC;
//			xMax_CD = xD;
//		} else {
//			xMin_CD = xD;
//			xMax_CD = xC;
//		}
//		if (yC < yD) {
//			yMin_CD = yC;
//			yMax_CD = yD;
//		} else {
//			yMin_CD = yD;
//			yMax_CD = yC;
//		}
//		if (xMax_AB < xMin_CD || xMin_AB > xMax_CD || yMax_AB < yMin_CD || yMin_AB > yMax_CD)
//			return null;
		
		if (xA < xB) {
			if (xC < xD) {
				if (xB<xC || xA>xD) return null;
			} else {
				if (xB<xD || xA>xC) return null;
			}
		} else {
			if (xC < xD) {
				if (xA<xC || xB>xD) return null;
			} else {
				if (xA<xD || xB>xC) return null;
			}
		}
		if (yA < yB) {
			if (yC < yD) {
				if (yB<yC || yA>yD) return null;
			} else {
				if (yB<yD || yA>yC) return null;
			}
		} else {
			if (yC < yD) {
				if (yA<yC || yB>yD) return null;
			} else {
				if (yA<yD || yB>yC) return null;
			}
		}
		

		double xI, yI;
		double K1, K2;

		if (xA==xB) {
			//------ 1er cas, le premier segment est vertical ------------------
			if (xC==xD)
				// les 2 segments verticaux
				return null;

			K2 = (yD-yC)/(xD-xC);
			xI = xA;
			yI = K2*(xI-xC) + yC;

			// xI,yI dans les segments [AB] et [CD] ?
			// TODO_KO cas trop rare pour nécessiter optimisation PERF : a optimiser
			if (yI >= Math.min(yA, yB) && yI <= Math.max(yA, yB) &&
					xI >= Math.min(xC, xD) && xI <= Math.max(xC, xD))
				return new Point_d(xI, yI, 5);
			else
				return null;

		} else if (xC==xD) {
			//------ 2eme cas, le deuxieme segment est vertical ----------------
			K1 = (yB-yA)/(xB-xA);
			xI = xC;
			yI = K1*(xI-xA) + yA;

			// xI,yI dans les segments [AB] et [CD] ?
			// TODO_KO cas trop rare pour nécessiter optimisation PERF : a optimiser les min max
			if (yI >= Math.min(yC, yD) && yI <= Math.max(yC, yD) &&
					xI >= Math.min(xA, xB) && xI <= Math.max(xA, xB))
				return new Point_d(xI, yI, 5);
			else
				return null;

		} else {
			//------ 3eme cas : cas standard -----------------------------------
			K1 = (yB-yA)/(xB-xA);
			K2 = (yD-yC)/(xD-xC);

			if (K1==K2) {
				// colinéarité -> 0 ou infinité de solution
				return null;
			}

			xI = (yC-K2*xC + K1*xA-yA) / (K1-K2);
			yI = K1*(xI-xA) + yA;

			// xI,yI dans les segments [AB] et [CD] ?
			// TODO_KO, pas reussi :-( PERF : a optimiser les min max
			if (xI >= Math.min(xA, xB) && xI <= Math.max(xA, xB) &&	xI >= Math.min(xC, xD) && xI <= Math.max(xC, xD))
				return new Point_d(xI, yI, 5);
//			if (xI >= xMin_AB && xI <= xMax_AB && xI >= xMin_CD && xI <= xMax_CD)
//				return new Point_d(xI, yI, 5);
			else
				return null;
		}
	}


	public static boolean doesPointBelongToSeg(Point_d pt, Segment_d seg, double toleranceSuperposition) {
		double minSegX, maxSegX, minSegY, maxSegY;
		if (seg.pt1.x < seg.pt2.x) {
			minSegX = seg.pt1.x;
			maxSegX = seg.pt2.x;
		} else {
			minSegX = seg.pt2.x;
			maxSegX = seg.pt1.x;
		}
		if (seg.pt1.y < seg.pt2.y) {
			minSegY = seg.pt1.y;
			maxSegY = seg.pt2.y;
		} else {
			minSegY = seg.pt2.y;
			maxSegY = seg.pt1.y;
		}
		
		// cas special
		if (seg.isVertical())
			return Math.abs(pt.x-seg.pt1.x) < toleranceSuperposition && pt.y>=minSegY && pt.y<=maxSegY;
		if (seg.isHorizontal())
			return Math.abs(pt.y-seg.pt1.y) < toleranceSuperposition && pt.x>=minSegX && pt.x<=maxSegX;
		
		//1er test simple : le pt doit etre dans le rect englobant du segment
		if (pt.x < minSegX || pt.x > maxSegX || pt.y < minSegY || pt.y > maxSegY)
			return false;
		
		// pour augmenter la precision des calculs et ne pas faire de division par des nombres trop petits
		//  on decoupe le probleme en 2
		// cela permet de comparer des pentes < 1 , necessaire a l'introduction d'un seuil
		if ((maxSegX-minSegX) > (maxSegY-minSegY)) {
			// le segment est plutot "horizontal"
			double dySeg_at_xPt = (maxSegY-minSegY) / (maxSegX-minSegX) * (pt.x - minSegX);
			double dyPt = pt.y - minSegY;
			return Math.abs(dySeg_at_xPt - dyPt) < toleranceSuperposition;
			
		} else {
			// le segment est plutot "vertical"
			double dxSeg_at_yPt = (maxSegX-minSegX) / (maxSegY-minSegY) * (pt.y - minSegY);
			double dxPt = pt.x - minSegX;
			return Math.abs(dxSeg_at_yPt - dxPt) < toleranceSuperposition;
		}
	}
	
	
	//**************************************************************************
	//**************************************************************************
	//**************************************************************************

	//Distance entre 2 rect
	static public double getDistRectRect (Rect_d r1, Rect_d r2) {
	    //CheckBoundaries r1
	    //CheckBoundaries r2

	    //on commence par calculer les distances suivant l'horizontale et la verticale
	    double dx;
	    double dy;

	    //Horizontalement
	    if (r2.xMin >= r1.xMax)
	        // r2 est à droite de r1
	        dx = r2.xMin - r1.xMax;
	    else if (r2.xMax <= r1.xMin)
	        // r2 est à gauche de r1
	        dx = r1.xMin - r2.xMax;
	    else
	        //sinon, r1 et r2 se superposent "horizontalement", la distance horizontale est 0
	        dx = 0;

	    //Verticalement
	    if (r2.yMin >= r1.yMax)
	        // r2 est au dessus de r1
	        dy = r2.yMin - r1.yMax;
	    else if (r2.yMax <= r1.yMin)
	        // r2 est au dessous de r1
	        dy = r1.yMin - r2.yMax;
	    else
	        //sinon, r1 et r2 se superposent "verticalement", la distance verticale est 0
	        dy = 0;

	    //D'ou la distance entre r1 et r2 (on fait qques tests pour éviter des appels à Sqr inutiles)
	    if (dx == 0)
	        return dy;
	    else if (dy == 0)
	        return dx;
	    else
	        return Math.sqrt (dx * dx + dy * dy);
	}

	// Détermine si 2 rectangles s'intersectent (en partie ou inclusion totale)
	// on utilise le meme principe que getDistRectRect() , mais sans calcul de distance
	static public boolean areRectIntersected (Rect_d r1, Rect_d r2) {
	    if (r2.xMin >= r1.xMax) return false; // r2 est à droite de r1
	    if (r2.xMax <= r1.xMin) return false; // r2 est à gauche de r1
	    if (r2.yMin >= r1.yMax) return false; // r2 est au dessus de r1
	    if (r2.yMax <= r1.yMin) return false; // r2 est au dessous de r1
	    //Arrivés la, on sait que r1 et r2 s'intersectent
	    return true;
	}

	// Surface d'intersection de 2 rectangles
	static public double getSurfaceRectIntersection (Rect_d r1, Rect_d r2) {
		if (areRectIntersected(r1, r2)) {
			double largIntersect = Math.min(r1.xMax, r2.xMax) - Math.max(r1.xMin, r2.xMin);
			double hautIntersect = Math.min(r1.yMax, r2.yMax) - Math.max(r1.yMin, r2.yMin);
//			assert(largIntersect>=0 && hautIntersect>=0); // obligatoire du fait du test positif d'intersection
			return largIntersect * hautIntersect;
		} else
			return 0;
	}

	//Distance entre un rect et un point
	static public double getDistRectPoint (Rect_d r, double x, double y) {
		//CheckBoundaries r

		//on commence par calculer les distances suivant l'horizontale et la verticale
	    double dx;
	    double dy;

	    //Horizontalement
	    if (x >= r.xMax)
	        // le point est à droite de r
	        dx = x - r.xMax;
	    else if (x <= r.xMin)
	        // le point est à gauche de r
	        dx = r.xMin - x;
	    else
	        //sinon, le point est "horizontalement" dans r, la distance horizontale est 0
	        dx = 0;

	    //Verticalement
	    if (y >= r.yMax)
	        // le point est au dessus de r
	        dy = y - r.yMax;
	    else if (y <= r.yMin)
	        // le point est en dessous de r
	        dy = r.yMin - y;
	    else
	        //sinon, le point est "verticalement" dans r, la distance verticale est 0
	        dy = 0;

	    //D'ou la distance entre le rect et le point (on fait qques tests pour éviter des appels à sqrt inutiles)
	    if (dx == 0)
	        return dy;
	    else if (dy == 0)
	    	return dx;
	    else
	    	return Math.sqrt (dx * dx + dy * dy);
	}

	// Distance min entre les coins d'un rect et un point
	static public double getDistCoinRectPoint (Rect_d r, double x, double y) {
		double dist_sq;
		//test des 4 coins
		dist_sq = (x-r.xMin)*(x-r.xMin) + (y-r.yMin)*(y-r.yMin);
		dist_sq = Math.min(dist_sq, (x-r.xMin)*(x-r.xMin) + (y-r.yMax)*(y-r.yMax));
		dist_sq = Math.min(dist_sq, (x-r.xMax)*(x-r.xMax) + (y-r.yMax)*(y-r.yMax));
		dist_sq = Math.min(dist_sq, (x-r.xMax)*(x-r.xMax) + (y-r.yMin)*(y-r.yMin));

		return Math.sqrt(dist_sq);
	}


	// Déplace le pt vers le plus proche emplacement situé sur le bord du rect
	static public void movePointToClosestBord_rect(Rect_d r, Point_d pt) {
		if (r.contientPt_strict(pt)) {
			int bordPlusProche;
			double distMin;// = Double.POSITIVE_INFINITY;
			double dist;// = 0;
			
			// 1) dist au bord gauche
			distMin = pt.x-r.xMin;
			bordPlusProche = 1;

			// 2) dist au bord haut
			dist = r.yMax - pt.y;
			if (dist < distMin) {
				distMin = dist;
				bordPlusProche = 2;
			}

			// 3) dist au bord droit
			dist = r.xMax - pt.x;
			if (dist < distMin) {
				distMin = dist;
				bordPlusProche = 3;
			}

			// 4) dist au bord droit
			dist = pt.y - r.yMin;
			if (dist < distMin) {
				//distMin = dist;
				bordPlusProche = 4;
			}

			// On prend le bord le plus proche
			if (bordPlusProche == 1) {
				pt.x = r.xMin;
			} else if (bordPlusProche == 2) {
				pt.y = r.yMax;
			} else if (bordPlusProche == 3) {
				pt.x = r.xMax;
			} else  /* if (bordPlusProche == 4) */ {
				pt.y = r.yMin;
			}

		} else {
			if (pt.x >= r.xMax)
				pt.x = r.xMax;
			else if (pt.x <= r.xMin)
				pt.x = r.xMin;
			// sinon, on laisse le x inchangé
			
			if (pt.y >= r.yMax)
				pt.y = r.yMax;
			else if (pt.y <= r.yMin)
				pt.y = r.yMin;
			// sinon, on laisse le y inchangé
		}
	}


	//**************************************************************************
	//**************************************************************************
	//**************************************************************************
	//static ArrayList<Integer> vIndexPointToKeep = new ArrayList<Integer>();
	public static ArrayList<Point_d> douglasPeuckerReduction_global (ArrayList<Point_d> vPoint, double tolerance) {

		if (vPoint.size() <= 2)
			return vPoint;

		int firstPoint = 0;
		int lastPoint = vPoint.size()-1;

		TIntArrayList vIndexPointToKeep = new TIntArrayList(vPoint.size());
		//vIndexPointToKeep.clear();

		// On ajoute le 1er et le dernier point
		vIndexPointToKeep.add (firstPoint);
		vIndexPointToKeep.add (lastPoint);

		douglasPeuckerReduction_recurs (vPoint, firstPoint, lastPoint, tolerance, vIndexPointToKeep);

		
		// Petite glut : si le 1er point = dernier point et que la reduction ne laisse plus que 3 points, alors
		//  il resterait 2 segments superposés, ce qui nous genera dans la suite du traitement
		if (vIndexPointToKeep.size()==3 && vPoint.get(0).equalsXY(vPoint.get(vPoint.size()-1))) {
			vIndexPointToKeep.clear();
			vIndexPointToKeep.add (firstPoint);
			vIndexPointToKeep.add (lastPoint);
			douglasPeuckerReduction_recurs (vPoint, firstPoint, lastPoint, tolerance/2, vIndexPointToKeep);
			if (vIndexPointToKeep.size()==3) {
				vIndexPointToKeep.clear();
				vIndexPointToKeep.add (firstPoint);
				vIndexPointToKeep.add (lastPoint);
				douglasPeuckerReduction_recurs (vPoint, firstPoint, lastPoint, tolerance/4, vIndexPointToKeep);
					if (vIndexPointToKeep.size()==3)
						return vPoint; // on abandonne et on retourne le vecteur de point inchangé ...
			}
		}
			
		
		// Attention l'appel précédent renvoi un vecteur d'index non trié !
		vIndexPointToKeep.sort();

		// Enfin, à partir des points conservés par l'algo, on reconstitue une liste de points, que l'on ajoute à vPointResult
		ArrayList<Point_d> vPointResult = new ArrayList<Point_d>(vIndexPointToKeep.size());
		Point_d pt;
		int index;
		int nbPointToKeep = vIndexPointToKeep.size();
		for (int i=0 ; i<nbPointToKeep ; i++) {
			index = vIndexPointToKeep.get(i);
			pt = vPoint.get(index);
			vPointResult.add (pt);
		}
		return vPointResult;
	}


	public static void douglasPeuckerReduction_recurs (ArrayList<Point_d> vPoint, int firstPoint, int lastPoint, double tolerance, TIntArrayList vIndexPointToKeep) {

		double distMax_sq = 0;
		int indexFarthest = 0;

		// Fin de la récursion
		if (firstPoint>=lastPoint-1)
			return;

		// Le 1er et le dernier point ne peuvent être identiques
		while (vPoint.get(firstPoint).equalsXY(vPoint.get(lastPoint))) {
			lastPoint--;
			vIndexPointToKeep.add (lastPoint);
			if (firstPoint>=lastPoint-1) return;
		}

		// Recherche de l'index du point le plus loin du segment [firstPoint, lastPoint]
		Point_d ptA, ptB, ptI;
		ptA = vPoint.get(firstPoint);
		ptB = vPoint.get(lastPoint);
		
		for (int i=firstPoint+1 ; i<lastPoint ; i++) {
			ptI = vPoint.get(i);

			double dist_sq = getDistancePtSegment_sq(ptA, ptB, ptI);

			if (dist_sq>distMax_sq) {
				distMax_sq = dist_sq;
				indexFarthest = i;
			}
		}

		// Si ce point le plus loin dépasse la tolérence, on le conserve et on appelle l'algo sur les 2 bouts de lignes autour.
		if (distMax_sq>(tolerance*tolerance) && indexFarthest!=0) {
			vIndexPointToKeep.add(indexFarthest);
			douglasPeuckerReduction_recurs (vPoint, firstPoint, indexFarthest, tolerance, vIndexPointToKeep);
			douglasPeuckerReduction_recurs (vPoint, indexFarthest, lastPoint, tolerance, vIndexPointToKeep);
		}
	}


	//**************************************************************************
	//**************************************************************************
	//**************************************************************************


	/**
	 * Contenance large (dans l'ellipse ou sur son pourtour
	 */
	public static boolean ellipseContientPoint_larg(Rect_d ellipse, Point_d pt) {
		// Paramètres de l'ellipse
		double a, b;
		a = ellipse.getWidth()/2;
		b = ellipse.getHeight()/2;
		Point_d ptCenterEllipse = ellipse.getCenter();

		return  ((pt.x-ptCenterEllipse.x)/a) * ((pt.x-ptCenterEllipse.x)/a) +
				((pt.y-ptCenterEllipse.y)/b) * ((pt.y-ptCenterEllipse.y)/b) <= 1;
	}
	public static boolean ellipseContientPoint_strict(Rect_d ellipse, Point_d pt) {
		// Paramètres de l'ellipse
		double a, b;
		a = ellipse.getWidth()/2;
		b = ellipse.getHeight()/2;
		Point_d ptCenterEllipse = ellipse.getCenter();

		return  ((pt.x-ptCenterEllipse.x)/a) * ((pt.x-ptCenterEllipse.x)/a) +
				((pt.y-ptCenterEllipse.y)/b) * ((pt.y-ptCenterEllipse.y)/b) < 1;
	}

	public static int getNbIntersect_Ellipse_Segment (Rect_d ellipse, Segment_d seg) {
		double[] tabCoord = new double[4];
		return computeIntersect_Ellipse_Segment(ellipse, seg, tabCoord);
	}

	/**
	 * Calcule les points d'intersection entre une ellipse et un segment et les stocke dans tabCoord
	 * Ne renvoie pas les intersections si elles sont extrémité du segment
	 * @param ellipse doit être non réduite à un segment, sinon erreur bloquante
	 * @param seg
	 * @param tabCoord doit contenir au moins la place pour 4 double
	 * @return le nombre de pt d'intersection
	 */
	public static int computeIntersect_Ellipse_Segment (Rect_d ellipse, Segment_d seg, double[] tabCoord) {
		assert ellipse.getWidth()>0 && ellipse.getHeight()>0;

		// Bounding box du segment
		double xSegMin, xSegMax, ySegMin, ySegMax;
		if (seg.pt1.x < seg.pt2.x) {
			xSegMin = seg.pt1.x;
			xSegMax = seg.pt2.x;
		} else {
			xSegMin = seg.pt2.x;
			xSegMax = seg.pt1.x;
		}
		if (seg.pt1.y < seg.pt2.y) {
			ySegMin = seg.pt1.y;
			ySegMax = seg.pt2.y;
		} else {
			ySegMin = seg.pt2.y;
			ySegMax = seg.pt1.y;
		}
		
		//**********************************************************************
		// 0) Cas triviaux
		//  a) segment hors du rect englobant de l'ellipse
		if (ySegMin > ellipse.yMax) return 0; //segment au dessus du rect
	    if (ySegMax < ellipse.yMin) return 0; //segment au dessous du rect
	    if (xSegMin > ellipse.xMax) return 0; //segment à droite du rect
	    if (xSegMax < ellipse.xMin) return 0; //segment à gauche du rect
		
		//  b) Segment réduit à un point : par convention, on ne renvoie rien (car le but final est de découper des segments ...)
		if (seg.getLengthSq()==0.0)
			return 0;


		//**********************************************************************
		// 1) Paramètres de l'ellipse
		double a, b;
		a = ellipse.getWidth()/2;
		b = ellipse.getHeight()/2;
		Point_d ptCenterEllipse = ellipse.getCenter();
		
		//**********************************************************************
		// 2) Paramètres de la droite : (u,v,w) /  u.x + v.y = w
		//   a) tout d'abord, u,v,w dans le repère d'origine
		double u,v,w;
		if (seg.isVertical()) {
			u = 1.0;
			v = 0.0;
			w = seg.pt1.x;
		} else if (seg.isHorizontal()) {
			u = 0.0;
			v = 1.0;
			w = seg.pt1.y;
		} else {
			v = 1.0; // convention
			w = (seg.pt1.y * seg.pt2.x - seg.pt2.y * seg.pt1.x) / (seg.pt2.x - seg.pt1.x);
			// x1 et x2 ne peuvent etre nuls simultanément, on choisit celui qui n'est pas nul
			if (seg.pt1.x != 0.0)
				u = (w - seg.pt1.y) / seg.pt1.x;
			else
				u = (w - seg.pt2.y) / seg.pt2.x;

		}
		//  b) puis translation dans le repère de l'ellipse
		w = w - u*ptCenterEllipse.x - v*ptCenterEllipse.y;

		//**********************************************************************
		// 3) Résolution avec intersection ligne
		int nbSol = computeIntersect_EllipseCentree_Line (a, b, u, v, w, tabCoord);
		if (nbSol >= 1) {
			tabCoord[0] += ptCenterEllipse.x; // on se remet dans le repère d'origine
			tabCoord[1] += ptCenterEllipse.y; // on se remet dans le repère d'origine
//			System.out.println(" x1 = " + tabCoord[0]);
//			System.out.println(" y1 = " + tabCoord[1]);
		}
		if (nbSol == 2) {
			tabCoord[2] += ptCenterEllipse.x; // on se remet dans le repère d'origine
			tabCoord[3] += ptCenterEllipse.y; // on se remet dans le repère d'origine
//			System.out.println(" x2 = " + tabCoord[2]);
//			System.out.println(" y2 = " + tabCoord[3]);
		}

		//**********************************************************************
		// 4) Restriction des solutions à l'intérieur du segment
		boolean bPt1OK = false;
		boolean bPt2OK = false;
		if (nbSol >= 1) {
			double x = tabCoord[0];
			double y = tabCoord[1];
			// on teste sur les x ou les y (pour les cas de segment horizontaux ou verticaux)
			if ( (x > seg.pt1.x && x < seg.pt2.x) || (x > seg.pt2.x && x < seg.pt1.x) ||
				 (y > seg.pt1.y && y < seg.pt2.y) || (y > seg.pt2.y && y < seg.pt1.y) ) {
				bPt1OK = true;
			} else {
				bPt1OK = false;
			}
		}
		if (nbSol == 2) {
			double x = tabCoord[2];
			double y = tabCoord[3];
			// on teste sur les x ou les y (pour les cas de segment horizontaux ou verticaux)
			if ( (x > seg.pt1.x && x < seg.pt2.x) || (x > seg.pt2.x && x < seg.pt1.x) ||
				 (y > seg.pt1.y && y < seg.pt2.y) || (y > seg.pt2.y && y < seg.pt1.y) ) {
				bPt2OK = true;
			} else {
				bPt2OK = false;
			}
		}

		if (bPt1OK && bPt2OK)
			return 2;
		else if (bPt1OK && !bPt2OK)
			return 1;
		else if (!bPt1OK && bPt2OK) {
			tabCoord[0] = tabCoord[2];
			tabCoord[1] = tabCoord[3];
			return 1;
		} else // (!bPt1OK && !bPt2OK)
			return 0;
	}


	private static int computeIntersect_EllipseCentree_Line (double a, double b, double u, double v, double w, double[] tabCoord) {
		double x, y;

		if (u != 0.0) {
			// équation du second degré
			double A, B, C;
			double temp = 1/(u*u*a*a);
			A = v*v*temp + 1/(b*b);
			B = -2*w*v*temp;
			C = w*w*temp - 1;

			assert !(A==0 && B==0 && C==0); // c'est censé être impossible

			if (A==0 && B==0 && C!=0)
				// pas de solution
				return 0;

			else if (A==0 && B!=0) {
				// 1 solution
				y = -C/B;
				x = (w-v*y)/u;
				tabCoord[0] = x;
				tabCoord[1] = y;
				return 1;

			} else if (A!=0) {
				double discrim = B*B - 4*A*C;
				//System.out.println("discrim = " + discrim);
				//System.out.println("A=" + A + "  B=" + B + " C=" + C);
				if (discrim>0) {
					// 2 solutions
					y = -(Math.sqrt(discrim) + B) / (2*A);
					x = (w-v*y)/u;
					tabCoord[0] = x;
					tabCoord[1] = y;
					y = (Math.sqrt(discrim) - B) / (2*A);
					x = (w-v*y)/u;
					tabCoord[2] = x;
					tabCoord[3] = y;
					return 2;
				} else if (discrim==0) {
					// 1 solution
					y = - B / (2*A);
					x = (w-v*y)/u;
					tabCoord[0] = x;
					tabCoord[1] = y;
					return 1;

				} else {
					// que des solutions complexes
					return 0;
				}


			} else {
				// impossible, on a déjà tout testé, normalement
				assert false;
				return 0;
			}

		} else {
			// Cas u=0 (segment horizontal)
			assert (v != 0);

			y = w/v;

			double temp = w*w/(v*v*b*b);
			if (1 == temp) {
				// 1 solution
				x = 0;
				tabCoord[0] = x;
				tabCoord[1] = y;
				return 1;
			} else if (1 > temp) {
				// 2 solutions
				x = a*Math.sqrt(1-temp);
				tabCoord[0] = x;
				tabCoord[1] = y;
				x = -x;
				tabCoord[2] = x;
				tabCoord[3] = y;
				return 2;
			} else {
				// pas de solution
				return 0;
			}
		}
	}

	// Déplace le pt vers le plus proche emplacement situé sur le bord de l'ellipse
	static public void movePointToClosestBord_ellipse(Rect_d r, Point_d pt) {

		Vector_d vect = new Vector_d(r.getCenter(), pt);

		if (vect.isPoint()) {
			// on est au centre de l'ellipse, on bouge le point de facon arbitraire sur le bord
			if (r.xMax-r.xMin > r.yMax-r.yMin)
				pt.y = r.yMax;
			else
				pt.x = r.xMax;
		}
		
		// le vecteur est agrandi pour etre au moins aussi lg que le plus grand cote de l'ellipse
		//  -> son extrémité sera forcément hors de l'ellipse
		double largestCoteEllipse = Math.max(r.xMax-r.xMin, r.yMax-r.yMin);

		vect.multiplyBy(largestCoteEllipse / vect.getLength());

		double[] tabCoord = new double[4];
		int nbIntersect = computeIntersect_Ellipse_Segment(r, vect, tabCoord);
		if (nbIntersect == 1) {
			pt.x = tabCoord[0];
			pt.y = tabCoord[1];
		} else
			assert (false);
	}


	//**************************************************************************
	//**************************************************************************
	//**************************************************************************

	// Détermine si 2 rectangles s'intersectent (en partie ou inclusion totale)
	// on utilise le meme principe que getDistRectRect() , mais sans calcul de distance
	static public boolean areRectSegIntersected (Rect_d r, Segment_d seg) {
		return computeSegIntersect_Rect_Segment(r,seg) != null;
	}
	static public boolean areRectSegIntersected__fast (Rect_d r, Segment_d seg, double[] tabCoord) {
		return computeLgIntersect_Rect_Segment__fast(r, seg, tabCoord) > 0;
	}

	static public double computeLgIntersect_Rect_Segment (Rect_d rect, Segment_d seg) {
		double[] tabCoord = new double[4];
		return computeLgIntersect_Rect_Segment__fast(rect, seg, tabCoord);
	}
	static public double computeLgIntersect_Rect_Segment__fast (Rect_d rect, Segment_d seg, double[] tabCoord) {

//		Segment_d segIntersect = computeSegIntersect_Rect_Segment(rect, seg);
//		if (segIntersect==null)
//			return 0;
//		else
//			return segIntersect.getLength();

		// MOINS PROPRE MAIS CONSOMME BCP MOINS DE MEMOIRE SI TRES NOMBREUX APPELS
		int nbPtIntersect = computeIntersect_Rect_Segment(rect, seg, tabCoord);

	    if (nbPtIntersect == 0) {
	        //si le seg est à l'intérieur du rect, alors l'intersection sera le segment lui-même
	        //sinon, l'intersection sera nulle
	        //if (rect.contientPt_larg(seg.getCenter())) // ==> TROP LENT
			if (rect.contientPt_larg(seg.pt1))
	            // intersection = [pt1,pt2]
	            return seg.getLength();
	        else
	        	return 0.0;

	    } else if (nbPtIntersect == 1) {
	        //il faut trouver le point du segment qui se trouve à l'intérieur du rect
	        if (rect.contientPt_strict(seg.pt1))
	            // intersection = [pt1,I]
				return Math.sqrt((seg.pt1.x-tabCoord[0])*(seg.pt1.x-tabCoord[0]) + (seg.pt1.y-tabCoord[1])*(seg.pt1.y-tabCoord[1]));
	        else if (rect.contientPt_strict(seg.pt2))
	            // intersection = [pt2,I]
				return Math.sqrt((seg.pt2.x-tabCoord[0])*(seg.pt2.x-tabCoord[0]) + (seg.pt2.y-tabCoord[1])*(seg.pt2.y-tabCoord[1]));
	        //sinon, c'est un cas particulier, ex. l'un des point a l'ext du rect, et l'autre sur un segment
	        else
	        	return 0.0;
	    } else {
	        //2 points d'intersection : intersection = [I,J];
	    	return Math.sqrt((tabCoord[0]-tabCoord[2])*(tabCoord[0]-tabCoord[2]) + (tabCoord[1]-tabCoord[3])*(tabCoord[1]-tabCoord[3]));
	    }
	}
	

	static public Segment_d computeSegIntersect_Rect_Segment (Rect_d rect, Segment_d seg) {
		double[] tabCoord = new double[4];
		Segment_d seg_out = new Segment_d();
		return computeSegIntersect_Rect_Segment__fast(rect, seg, tabCoord, seg_out);
	}
	// tabCoord constitue une zone memoire pre-allouée ou seront ecrit les coord de l'intersection
	// seg_out est un segment pre-aloué ou sera stocké les coordonées du segment intersection
	// il faut regarder le segment retourné par la fonction ! qui peut etre different de seg_out
	// car si pas d'intersection, le retour est null et seg_out est inchangé !!
	static public Segment_d computeSegIntersect_Rect_Segment__fast (Rect_d rect, Segment_d seg_in, double[] tabCoord, Segment_d seg_out) {
		
		int nbPtIntersect = computeIntersect_Rect_Segment(rect, seg_in, tabCoord);

	    if (nbPtIntersect == 0) {
	        //si le seg est à l'intérieur du rect, alors l'intersection sera le segment lui-même
	        //sinon, l'intersection sera nulle
	        //if (rect.contientPt_larg(seg.getCenter())) // ==> TROP LENT
			if (rect.contientPt_larg(seg_in.pt1))
	            // intersection = [pt1,pt2]
	            return seg_in;
	        else
	        	return null;

	    } else if (nbPtIntersect == 1) {

			//il faut trouver le point du segment qui se trouve à l'intérieur du rect
	        if (rect.contientPt_strict(seg_in.pt1)) {
	            // intersection = [pt1,I]
				//return new Segment_d(seg_in.pt1, new Point_d(tabCoord[0], tabCoord[1]));
				seg_out.setCoord(seg_in.pt1.x, seg_in.pt1.y, tabCoord[0], tabCoord[1]);
				return seg_out;
	        	
			} else if (rect.contientPt_strict(seg_in.pt2)) {
	            // intersection = [pt2,I]
	        	//return new Segment_d(seg_in.pt2, new Point_d(tabCoord[0], tabCoord[1]));
				seg_out.setCoord(seg_in.pt2.x, seg_in.pt2.y, tabCoord[0], tabCoord[1]);
				return seg_out;

			//sinon, c'est un cas particulier, ex. l'un des point a l'ext du rect, et l'autre sur un segment
			} else
	        	return null;
	    } else {
	        //2 points d'intersection : intersection = [I,J];
	    	//return new Segment_d(new Point_d(tabCoord[0], tabCoord[1]), new Point_d(tabCoord[2], tabCoord[3]));
			seg_out.setCoord(tabCoord[0], tabCoord[1], tabCoord[2], tabCoord[3]);
			return seg_out;
	    }
	}


	/**
	 * @param rect NE DOIT pas etre reduit a un point ou segment (sinon erreur bloquante)
	 * @param seg
	 * @param tabCoord
	 * @return on ne renvoie rien si intersection au niveau de l'extrémité du segment, et rien si rect réduit à segment
	 */
	static public int computeIntersect_Rect_Segment (Rect_d rect, Segment_d seg, double[] tabCoord) {

		assert (rect.xMin!=rect.xMax && rect.yMin!=rect.yMax);

		boolean bSegVertical = seg.isVertical();
		boolean bSegHorizontal = seg.isHorizontal();

		// Bounding box du segment
		// TODO_PERF_ANNULE_CAR_PEU_LISIBLE : faire les tests de hors rect directement ici pour gagner plus remonter tout ce code des le debut
		double xSegMin, xSegMax, ySegMin, ySegMax;
		if (seg.pt1.x < seg.pt2.x) {
			xSegMin = seg.pt1.x;
			xSegMax = seg.pt2.x;
		} else {
			xSegMin = seg.pt2.x;
			xSegMax = seg.pt1.x;
		}
		if (seg.pt1.y < seg.pt2.y) {
			ySegMin = seg.pt1.y;
			ySegMax = seg.pt2.y;
		} else {
			ySegMin = seg.pt2.y;
			ySegMax = seg.pt1.y;
		}

		//**********************************************************************
		// 0) Cas triviaux
		//  a) Segment réduit à un point : par convention, on ne renvoie rien (car le but final est de découper des segments ou de calculer des longeur d'intersection)
		if (bSegVertical && bSegHorizontal)	return 0;
	    //  b) on élimine dès le début le cas fréquent ou le segment est completement hors du rect
	    if (ySegMin >= rect.yMax) return 0; //segment au dessus du rect
	    if (ySegMax <= rect.yMin) return 0; //segment au dessous du rect
	    if (xSegMin >= rect.xMax) return 0; //segment à droite du rect
	    if (xSegMax <= rect.xMin) return 0; //segment à gauche du rect

	    int nbPtIntersect = 0;

	    double xI;
	    double yI;

		// Intersection avec le coté supérieur du rect
	    if (ySegMax > rect.yMax) { //sinon, pas d'intersection possible
			// Si segment horizontal -> pas d'intersection, car on ne compte pas la frontière du rect
	        if (!bSegHorizontal) {
	            //segment vertical ou oblique, on utilise notre cher Thalès
	            xI = (seg.pt1.y - rect.yMax) * (seg.pt2.x - seg.pt1.x) / (seg.pt1.y - seg.pt2.y) + seg.pt1.x;
	            yI = rect.yMax;

	            //on teste si ce pt d'intersection se trouve bien a l'interieur du coté du rect
	            // cf. la présence stricte du pt dans le segment est assurée par les tests précédents
	            if (xI >= rect.xMin && xI <= rect.xMax) {
					tabCoord[0] = xI;
	                tabCoord[1] = yI;
					nbPtIntersect++;
	        	}
	        }
		}

	    //intersection avec le coté inférieur du rect
	    if (ySegMin < rect.yMin) { //sinon, pas d'intersection possible
	        // Si segment horizontal -> pas d'intersection, car on ne compte pas la frontière du rect
	        if (!bSegHorizontal) {
	            //segment vertical ou oblique, on utilise notre cher Thalès
	            xI = (seg.pt1.y - rect.yMin) * (seg.pt2.x - seg.pt1.x) / (seg.pt1.y - seg.pt2.y) + seg.pt1.x;
	            yI = rect.yMin;

	            //on teste si ce pt d'intersection se trouve bien a l'interieur du coté du rect
	            // cf. la présence stricte du pt dans le segment est assurée par les tests précédents
	            if (xI >= rect.xMin && xI <= rect.xMax) {
					if (nbPtIntersect == 0) {
						tabCoord[0] = xI;
						tabCoord[1] = yI;
					} else { // nbPtIntersect == 1
						tabCoord[2] = xI;
						tabCoord[3] = yI;
					}
					nbPtIntersect++;
	            }
	        }
	    }


	    //intersection avec le coté gauche du rect
	    if (nbPtIntersect < 2) {
	        if (xSegMin < rect.xMin) { //sinon, pas d'intersection possible
	            // Si segment vertical -> pas d'intersection, car on ne compte pas la frontière du rect
				if (!bSegVertical) {
	                //segment horizontal ou oblique, on utilise notre cher Thalès
	                xI = rect.xMin;
	                yI = (seg.pt1.x - rect.xMin) * (seg.pt2.y - seg.pt1.y) / (seg.pt1.x - seg.pt2.x) + seg.pt1.y;

					//on teste si ce pt d'intersection se trouve bien a l'interieur du coté du rect
					// cf. la présence stricte du pt dans le segment est assurée par les tests précédents
	                if (yI >= rect.yMin && yI <= rect.yMax) {

						// on teste si le point n'a pas déjà été trouvé sur le cote haut ou bas
						if (nbPtIntersect == 1 && tabCoord[0] == xI && tabCoord[1] == yI)
							;//System.out.println("pt deja trouvé");
						else {
							if (nbPtIntersect == 0) {
								tabCoord[0] = xI;
								tabCoord[1] = yI;
							} else { // nbPtIntersect == 1
								tabCoord[2] = xI;
								tabCoord[3] = yI;
							}
							nbPtIntersect++;
						}
	                }
	            }
	        }
	    }

	    //intersection avec le coté droit du rect
	    if (nbPtIntersect < 2) {
	        if (xSegMax > rect.xMax) { //sinon, pas d'intersection possible
				// Si segment vertical -> pas d'intersection, car on ne compte pas la frontière du rect
				if (!bSegVertical) {
	                //segment horizontal ou oblique, on utilise notre cher Thalès
	                xI = rect.xMax;
	                yI = (seg.pt1.x - rect.xMax) * (seg.pt2.y - seg.pt1.y) / (seg.pt1.x - seg.pt2.x) + seg.pt1.y;

	                //on teste si ce pt d'intersection se trouve bien a l'interieur du coté du rect
					// cf. la présence stricte du pt dans le segment est assurée par les tests précédents
	                if (yI >= rect.yMin && yI <= rect.yMax) {

						// on teste si le point n'a pas déjà été trouvé sur le cote haut ou bas
						if (nbPtIntersect == 1 && tabCoord[0] == xI && tabCoord[1] == yI)
							;//System.out.println("pt deja trouvé");
						else {
							if (nbPtIntersect == 0) {
								tabCoord[0] = xI;
								tabCoord[1] = yI;
							} else { // nbPtIntersect == 1
								tabCoord[2] = xI;
								tabCoord[3] = yI;
							}
							nbPtIntersect++;
						}
	                }
	            }
	        }
	    }
		return nbPtIntersect;
	}


//	public static void testAll_Intersect_Rect_Segment() {
//		System.out.println("---- TEST INTERSECT RECT / SEGMENT ----");
//
//		Segment_d seg;
//		Rect_d rect = new Rect_d(2,1,6,4);
//		// 1) Cas triviaux
//		testOne_Intersect_Rect_Segment ("special_0 haut", rect, new Segment_d(3,5,4,6), 0);
//		testOne_Intersect_Rect_Segment ("special_0 droite", rect, new Segment_d(7,2,12,3), 0);
//		testOne_Intersect_Rect_Segment ("special_0 bas", rect, new Segment_d(3,0,4,1), 0);
//		testOne_Intersect_Rect_Segment ("special_0 gauche", rect, new Segment_d(0,2.3,2,3.7), 0);
//		testOne_Intersect_Rect_Segment ("special_0 horizont", rect, new Segment_d(5,4,7,4), 0);
//		testOne_Intersect_Rect_Segment ("special_0 horizont", rect, new Segment_d(6,0,6,2), 0);
//
//		testOne_Intersect_Rect_Segment ("std_0 HG", rect, new Segment_d(2,4,3,3), 0);
//		testOne_Intersect_Rect_Segment ("std_0 HD", rect, new Segment_d(6,4,5,3), 0);
//		testOne_Intersect_Rect_Segment ("std_0 BDdiag", rect, new Segment_d(6,1,2,4), 0);
//		testOne_Intersect_Rect_Segment ("std_0 BG", rect, new Segment_d(2,1,3,2), 0);
//		testOne_Intersect_Rect_Segment ("std_0 ext", rect, new Segment_d(5,5,8,3), 0);
//		testOne_Intersect_Rect_Segment ("std_0 int", rect, new Segment_d(3,3,5,3.5), 0);
//
//		testOne_Intersect_Rect_Segment ("std_1 G hor", rect, new Segment_d(1,3,3,3), 1);
//		testOne_Intersect_Rect_Segment ("std_1 G obl", rect, new Segment_d(1,3,3,2), 1);
//		testOne_Intersect_Rect_Segment ("std_1 G obl extr", rect, new Segment_d(1,3,6,1), 1);
//		testOne_Intersect_Rect_Segment ("std_1 H vert", rect, new Segment_d(4,5,4,3), 1);
//		testOne_Intersect_Rect_Segment ("std_1 H obl", rect, new Segment_d(4,5,5,3), 1);
//		testOne_Intersect_Rect_Segment ("std_1 H obl extr", rect, new Segment_d(4,5,6,1), 1);
//		testOne_Intersect_Rect_Segment ("std_1 D hor", rect, new Segment_d(5,3,7,3), 1);
//		testOne_Intersect_Rect_Segment ("std_1 D obl", rect, new Segment_d(5,2,7,3), 1);
//		testOne_Intersect_Rect_Segment ("std_1 D obl extr", rect, new Segment_d(2,1,7,3), 1);
//		testOne_Intersect_Rect_Segment ("std_1 B vert", rect, new Segment_d(4,0,4,2), 1);
//		testOne_Intersect_Rect_Segment ("std_1 B obl", rect, new Segment_d(4,0,3,2), 1);
//		testOne_Intersect_Rect_Segment ("std_1 B obl extr", rect, new Segment_d(4,0,2,4), 1);
//
//		testOne_Intersect_Rect_Segment ("special_1 diag", rect, new Segment_d(6,1,-2,7), 1);
//		testOne_Intersect_Rect_Segment ("special_1 HD", rect, new Segment_d(7,5,5,3), 1);
//		testOne_Intersect_Rect_Segment ("special_1 BD", rect, new Segment_d(7,2,5,0), 1);
//		testOne_Intersect_Rect_Segment ("special_1 BG", rect, new Segment_d(1,0,3,2), 1);
//
//		testOne_Intersect_Rect_Segment ("std_2 hor", rect, new Segment_d(1,3,7,3), 2);
//		testOne_Intersect_Rect_Segment ("std_2 vert", rect, new Segment_d(3,0,3,5), 2);
//		testOne_Intersect_Rect_Segment ("std_2 diag1", rect, new Segment_d(-2,7,10,-2), 2);
//		testOne_Intersect_Rect_Segment ("std_2 diag2", rect, new Segment_d(-2,-2,10,7), 2);
//		testOne_Intersect_Rect_Segment ("std_2 corn BD", rect, new Segment_d(4,0,7,3), 2);
////		testOne_Intersect_Rect_Segment ("std_2 ", rect, new Segment_d(,,,), 2);
////		testOne_Intersect_Rect_Segment ("std_2 ", rect, new Segment_d(,,,), 2);
//	}
//
//	private static void testOne_Intersect_Rect_Segment(String testName, Rect_d rect, Segment_d seg, int nbIntersectExpected) {
//		double[] tabCoord = new double[4];
//		int result = UtilsGeom.computeIntersect_Rect_Segment(rect, seg, tabCoord);
//		if (result == nbIntersectExpected)
//			System.out.println( (testName + "                  ").substring(0, 18) + " : OK");
//		else
//			System.out.println( (testName + "                  ").substring(0, 18) + " : KO : result= "+result+", expected="+nbIntersectExpected);
//
//		if (result>=1)
//			System.out.println("x = " + tabCoord[0] + " , y = " + tabCoord[1]);
//		if (result==2)
//			System.out.println("x = " + tabCoord[2] + " , y = " + tabCoord[3]);
//	}

}
